A dynamical phase transition in a model for evolution with migration 
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Migration between different habitats is ubiquitous among biological populations. In this Letter, 
we study a simple quasispecies model for evolution in two different habitats, with different fitness 
landscapes, coupled through one-way migration. Our model applies to asexual, rapidly evolving 
organisms such as microbes. Our key finding is a dynamical phase transition at a critical value of 
the migration rate. The time to reach steady state diverges at this critical migration rate. Above 
the transition, the population is dominated by immigrants from the primary habitat. Below the 
transition, the genetic composition of the population is highly non-trivial, with multiple coexisting 
quasispecies which are not native to either habitat. Using results from localization theory, we show 
that the critical migration rate may be very small - demonstrating that evolutionary outcomes can 
be very sensitive to even a small amount of migration. 
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Biological dispersal — the movement of organisms be- 
tween habitats — is a ubiquitous phenomenon with im- 
portant and wide-ranging consequences. In the natu- 
ral environment, organisms expand their ranges, colonise 
new habitats, and can undergo speciation if they become 
spatially isolated. Dispersal plays a key role in deter- 
mining spatial and temporal patterns of genetic diver- 
sity in all organisms For sexual organisms, with 
low mutation rates, population subdivision into demes, 
connected by migration, can have important effects on 
genetic diversity [2, [H , while in continuous space, trans- 
mission of unfit alleles can prevent the expansion of a 
species' range For asexual, rapidly evolving organ- 
isms such as bacteria and viruses, dispersal also facilitates 
the emergence of new diseases and resistance to known 
treatments. The "source-sink" paradigm 0, @|, in which 
migration from a favourable habitat maintains organisms 
in an unfavourable one, has recently been used to ex- 
plain the microbial genetics of urinary tract infections Q • 
However, despite its importance, a general understand- 
ing of how migration affects mutation-selection balance 
in microbial systems is lacking. In particular, one would 
like to know how migration changes the proportions of 
different genotypes in the evolving population. 

In order to study the role of migration we introduce in 
this letter a simple statistical physics model for the evo- 
lutionary dynamics of migrating asexual organisms. Our 
model comprises two environmental habitats (with dif- 
ferent fitness landscapes) coupled by one-way migration 
of organisms from the primary to the secondary habitat. 
Using a quasispecies approach, we find that the model un- 
dergoes a dynamical phase transition: at a critical value 
of the migration rate, the time to reach the steady state 
diverges. For sub-critical migration rates, the steady- 
state population in the secondary habitat is made up of 
the organisms "native" (best adapted) to this habitat, as 



well as other, non-trivial, quasispecies, which are not na- 
tive to either habitat. Above the critical migration rate, 
the native quasispecies in the secondary habitat is wiped 
out by immigrants from the primary habitat. We use 
results from localization theory to gain insight into the 
transition and to show that the critical migration rate is 
typically small, demonstrating that even a small amount 
of migration can have an important effect on evolutionary 
dynamics. 

In our model, organisms have M possible genotypes. 
Ni and n,; denote the abundance (number density) of or- 
ganisms with genotype i in the primary and secondary 
habitat, respectively. The populations in the two habi- 
tats are thus described by the vectors N = (N±, . . . , Nm) 
and n = (m, . . . , Um)- Organisms migrate from the pri- 
mary to the secondary habitat with rate k. Within each 
habitat, mutations transform organism i to j with rate 
7-Aij, where Ay is a symmetric adjacency matrix, to be 
discussed later. Organisms of type i reproduce at a rate 
$j — Nj in the primary habitat and <j>i — ■ rij in the 

secondary habitat. The vectors $ = ($1, . . . , $a/) and 
4> = {4>x, . . . , 4>m) thus describe the fitness landscapes (or 
the maximal growth rate for organisms with genotype i) 
in each habitat. The terms — . Nj and — rij in the 
growth rates account for population saturation due to fi- 
nite resources, as in the logistic equation. This model is 
based on the para-mu-se (parallel mutation and selection) 
Q version of quasispecies theory widely discussed in 
the biological, chemical and physical literature [Iol - [l5| . 

The time evolution of the system is governed by the 
following set of equations for i = 1, . . . , M: 

AT = - n j) + 7 E A M ~ Ni), (1) 

o j 

hi = nj((f>i — rij) + 7 Ajj (rij - m) + kNj, (2) 

j j 
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where we have assumed that the primary habitat is large, 
so that the loss of individuals due to migration has a neg- 
ligible effect on its population [l6[. For the calculations 
presented here, we suppose that the fitness values 
are independent random numbers drawn from a distribu- 
tion P((p), common to both environments. Thus genomes 
which are well-adapted in the primary habitat are likely 
to be maladapted in the secondary habitat. 

We first present the analytical solution for the steady 
state . For the primary habitat it is known from qua- 
sispecies theory 0, 0| that the steady-state abundances 
N* are 



N* = 



Ai 
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with e = (1 
the matrix W t 



, , 1). Here Ax is the largest eigenvalue of 
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being the graph Laplacian), and 'J/i is the corresponding 
eigenvector. We now determine the steady-state geno- 
type abundances in the secondary habitat by expanding 
in the eigenbasis of Vu = + "fAij 



M f T jCt* 
^ That - A Q 
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where ip a and A Q are the eigenvectors and eigenvalues 
of Vtj (ordered as Ai > A2 > . . .) and n to t is the total 
steady-state population in the secondary habitat, which 
is determined self-consistently as the largest root of 
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To proceed further, we now make some specific as- 
sumptions about the structure of the genome space (the 
mutation matrix Aij) and the fitness landscape P(tp). To 
this end, we suppose that the mutation graph is a one di- 
mensional closed chain, in which mutations are possible 
only between neighbouring genotypes (i.e. Aij = 1 if 
i = (AI + j ± l)modM, and zero otherwise). We further 
suppose that the fitness can take only two values: 1 and 
with probability p and 1 — p, respectively. Since it has 
been suggested that viable genotypes form an intercon- 
nected network in genome space [18J, we shall consider 
the case p ~ 1, so that the fitness landscape is charac- 
terised by "islands" of fit genotypes separated by unfit 
ones. 

Figure Q] shows how the population composition in 
the secondary habitat depends on the migration rate k. 
When k is very small, the steady-state distribution n* is 
peaked around the longest sequence of maximal fitness 
values: this peak corresponds to the "native" (or best- 
adapted) quasispecies for the secondary habitat. When 
k is very large (much larger than the mutation rate 7), 
the secondary habitat becomes dominated by immigrants 
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Figure 1: Examples of steady-state genotype abundances in 
the secondary habitat n*, for different values of the migra- 
tion rate k where p = 0.9, 7 = 0.01 and M = 128. These 
results were obtained by numerical self-consistent solution of 
Eqs. ©, Q and The top left panel shows zero abun- 
dances (n* — 0) in the absence of migration. The inset shows 
the abundances N* in the primary habitat (red line). The 
top right panel shows the "native" steady-state n* (blue line), 
for very small migration rate, as well as the fitness landscape 
4>i/ 100 (black line). The other panels show n* in the sec- 
ondary habitat, for various values of the migration rate k. 
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Figure 2: Numerical results (from Eq. ([2]), using Eq. (0 for 
N*)) for the time T to reach the steady-state, starting from 
Ni = N* and n; = 0. Left panel: T as a function of migration 
rate k, for the same system as in Fig.[T] The steady state was 
assumed to have been reached when ^T \m(t + l)/rii(t) — 
1\/M < 10 -10 . Right panel: T(k/k ), where k is determined 
from Eq. @, for M = 64, 7 = 0.01, p = 0.7, normalized so 
that T(10) = 1. Results for 20 representative sets of <f> axe 
presented on a log-log plot. 



from the primary habitat and the distribution of n* tends 
to the primary-habitat steady-state distribution TV* . 

In contrast, for intermediate migration rates, the ge- 
netic composition in the secondary habitat is highly non- 
trivial. As k increases from zero, the quasispecies native 
to the secondary habitat is joined by additional, non- 
native, quasispecies peaks. These do not correspond to 
the native quasispecies from the primary habitat, but 
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Figure 3: Genetic composition in the secondary habitat rn(t), 
during the approach to steady state, for the same system as 
in Figure[T] for two migration rates k = 0.0001 (k <C ko, black 
line, left) and k — 0.003 (k S> ko, red, line right, the same 
vertical scale), for t — 1, . . . , 2 18 . 



are instead determined by the overlap of eigenvectors in 
the primary and secondary habitats (as in Eq. ((4])). As 
the migration rate is increased slightly further from 0.002 
to 0.003, these new peaks dominate completely and the 
native quasispecies of the secondary habitat disappears. 
This effect can be triggered by a very moderate change in 
the migration rate. The appearance of these new quasis- 
pecies peaks suggests that migration coupled to mutation 
can provide a mechanism for generation and maintenance 
of genetic diversity (as will be shown later in Figure @| . 
Finally, as the migration rate increases further, the new 
peaks merge into the native quasispecies peak from the 
primary habitat. 

Figure [2] (left panel) shows that this non-trivial depen- 
dence of the steady state population on the migration 
rate is accompanied by striking changes in the system 
dynamics. The time to reach the steady state plotted as 
a function of k shows a striking maximum at ko ~ 0.0027, 
suggesting a critical slowing down and a likely dynami- 
cal phase transition. The approach to the steady state 
for k <C ko is much slower than for k ^> fco . Figure [3] 
illustrates the underlying reason for this. Here we plot 
snapshots of n at various moments in time during the 
approach to the steady state, for the same parameter 
set and fitness landscape, for migration rates below and 
above fco. For both migration rates, the immigrating pop- 
ulation initially has the same composition as the primary 




habitat. For fc <C fco, the primary habitat quasispecies 
peak is lost entirely and the system undergoes a slow 
process of jumps between various local fitness maxima 
before finally settling in the global optimum. In con- 
trast, for k 3> fco; the system rapidly relaxes to a steady 
state which overlaps strongly with that of the primary 
habitat. 

Returning to our analytical expressions for n*, Eqs. (j4|) 
and ([5]), we can estimate the critical migration rate fco at 
which the dynamical phase transition takes place. Equa- 
tion ([3]) expresses n as a sum of eigenvectors ip for the 
secondary habitat, weighted by their overlap with N*. 
When fc — > 0, n to t — > M — 1 and ft — > tpi (see below). 
This is the native quasispecies solution for the secondary 
habitat [ljjj]. The phase transition occurs when this so- 
lution becomes dominated by the contributions from the 
other terms (a — 2, . . . ,M in Eqs. ((4]) and ([5])), which 
arise from overlap with the primary habitat solutions. 
This happens at a migration rate approximately given 
by 



(6) 



To show that this result indeed corresponds to the critical 
migration rate at which the transition happens, we plot 
in Figure [21 right panel, the time T to reach steady state 
as a function of fc/fco, where fco is determined from ([5]), 
for simulated dynamics on rj 20 representative random 
fitness sequences 0. Each fitness landscape generates a 
slightly different curve T(fc/fc ), but all the curves appear 
to diverge at fc = fc , indicating a phase transition. 

We now briefly discuss how the steady-state properties 
of the system are affected by this phase transition. For 
fc < fco, the a = 1 terms in Eqs. Q and §5§ dominate, 
while for fc > fco, the terms a = 2,...,M dominate. 
This has important consequences for the total population 
ntot in the secondary habitat: for fc < fco, ntot — 1 fl9| . 
while for fc > fco, n to t grows with fc. This prediction is 
confirmed numerically in Fig. 01 left panel, for a number 
of randomly generated fitness landscapes. 

Another important steady-state property is genetic di- 
versity. We noted from Figure [T] that multiple quasis- 
pecies peaks can coexist in steady state for intermediate 
migration rates. Figure [5J right panel, plots the partici- 
pation ratio (PR) r = n* ) 2 / X)i( n i ) 2 ■ as a function 
of fc/fco, for several realisations of the fitness landscape. 
The PR is a convenient measure of diversity, which shows 
how many of the n*'s are much larger than zero. Fig- 
ure 2] shows that r reaches a maximal value at about 
fc w 0.25 — 0.75fco and then decreases as fc approaches fco; 
above fco, the diversity remains approximately constant. 
Thus for weak migration the genetic diversity is increased 
whereas for strong migration it is washed out. 

For the one-dimensional model considered here, re- 
sults from localisation theory [2(| allow us to estimate 
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Figure 4: Left: Total steady-state population ntot in the sec- 
ondary habitat, as a function of rescaled migration rate k/ko, 
for M = 128, p = 0.9, 7 = 0.1, for a number of typical random 
fitness landscapes. Right: The normalized participation ratio 
r/ro for ro = r(k — > 0), for the same set of fitness landscapes. 



the value of the critical migration rate ko- The eigen- 
vector equation for "Jq in the primary habitat maps 
onto a Schrodinger equation with random potential Uj = 



(A^ a )j + Uj(^ a )j = E a (^ a )j 



(7) 



where E a = — A Q /7, and likewise for an eigenvector 
tp a in the secondary habitat. Equation f7|) is essen- 
tially a ID tight-binding electron model [2l|, in which 
Uj = — 1/7 with probability p and Uj = with prob- 
ability 1 — p. Localization theory tells us that for this 
problem the ground state eigenvector is localized, tak- 
ing the form ^ij ~ sin(j7r/?ii) on the longest run w of 



consecutive sites with Uj 



-I/7, and has eigenvalue 



Ai ~ 1 — 77r 2 /w 2 . Eigenvectors corresponding to ex- 
cited states are similarly localized on other, shorter po- 
tential wells. To estimate ko, we observe that the largest 
contribution to the sum in ^ comes from the eigenvec- 
tor with the greatest overlap with TV*, which we denote 
ipp. Assuming that N* and ipp are localized on poten- 
tial wells of length w and v, respectively, we can esti- 
mate that (ipj ■ N*)(ipJ ■ e) ~ v/w. The lengths w,v 
are the longest runs of Uj = — 1/7 in sequences of in- 
dependent binary random numbers of length M and w, 
respectively, therefore w ~ ln(M(l — p))/\n(l/p) and 
v ~ ln(iu(l —p))/ ln(l/p). For large M, v is much smaller 
than w, so Ai — \p ~ jtt 2 /v 2 . Inserting this into Eq. ©, 
and setting e = 1 — p, we finally obtain 



(fc ) ~ jw/v 



3 ^ 7 e 2 ln(Afe) 
(lnlnMe) 3 



(8) 



Remarkably, this rough estimate agrees up to a factor 
w 2 with our simulation results. Here we have considered 
small e, where multiple fit genomes lie close together in 
genotype space, and we see from ([8]) that (ko) is much 
smaller than 7 for moderately large Me. This means that 
even a very small migration rate (smaller than the muta- 
tion rate) can dramatically change the course of evolution 
in the secondary environment [22j. 



In summary, we have introduced a simple model for 
the evolution of asexual organisms in two coupled habi- 
tats with different fitness landscapes. We have shown 
that a dynamical phase transition occurs as the migra- 
tion rate changes. Bifurcations caused by migration have 
been observed in several models of sexual populations 
0il3| but, to our knowledge, the present work is the first 
to consider the effects of migration on the evolutionary 
dynamics of asexual organisms from a quasispecies per- 
spective. In our model, at the critical migration rate, the 
population in the secondary habitat becomes dominated 
by immigrants from the primary habitat. For subcritical 
migration rates, our quasispecies model also reveals that 
migration can provide a novel mechanism for creation 
and maintenance of genetic diversity. 

To obtain analytical results and clear insights into the 
physics of the model, we have mainly considered a simple 
one-dimensional closed-chain representation of the geno- 
type space and binary random fitness landscapes. As a 
step towards more complex and realistic representations 
of the genome space and fitness landscape, we have also 
carried out numerical simulations for a continuous, uni- 
form distribution of the fitness, as well as a hypercubic 
mutation graph. Our key results (in particular the dy- 
namical phase transition as a function of migration rate) 
remain valid in these cases, suggesting that our findings 
are likely to be of general significance. It will be inter- 
esting to extend our work to empirical fitness landscapes 
generated from experimental data (2^ .and, inspired by 
existing models for sexual organisms 0, 0|, and recent 
models in microbial ecology [24]], to multiple connected 
habitats and spatially varying environments. 
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